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We analyze the ground state phase diagram of attractive lattice bosons, which are stabilized by a 
three-body onsite hardcore constraint. A salient feature of this model is an Ising type transition from 
a conventional atomic superfluid to a dimer superfluid with vanishing atomic condensate. The study 
builds on an exact mapping of the constrained model to a theory of coupled bosons with polynomial 
interactions, proposed in a related paper [11]. In this framework, we focus by analytical means 
on aspects of the phase diagram which are intimately connected to interactions, and are thus not 
accessible in a mean field plus spin wave approach. First, we determine shifts in the mean field phase 
border, which are most pronounced in the low density regime. Second, the investigation of the strong 
coupling limit reveals the existence of a "continuous supersolid" , which emerges as a consequence of 
enhanced symmetries in this regime. We discuss its experimental signatures. Third, we show that 
the Ising type phase transition, driven first order via the competition of long wavelength modes at 
generic fillings, terminates into a true Ising quantum critical point in the vicinity of half filling. 
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I. INTRODUCTION 

It was recently recognised that two-body and three- 
body loss processes for bosons in an optical lattice could 
give rise to effective models involving two-body and 
three-body hardcore constraints, respectively. The two- 
body case was observed in an experiment with Feshbach 
molecules [TJ [5] , while it has been proposed theoretically 
to take advantage of strong three-body loss to create 
a three-body hardcore constraint in bosonic [3l 0] and 
fermionic [5] lattice systems. The mechanism behind the 
constraint is that the dissipative process suppresses co- 
herent tunnelling processes that would create double or 
triple occupation and lead to loss. 

A salient feature of a bosonic lattice gas with three- 
body onsite constraint is the possibility to tunc it to 
attractive two-body interactions. The associated dimer 
bound state formation has a profound effect on the many- 
body system, resulting in an Ising-type quantum phase 
transition from a conventional atomic superfluid to a 
dimer superfluid with vanishing atomic order parame- 
ter but nonzero pairing correlation. The possibility of 
observing Ising type behavior in cold atomic gases has 
been uncovered earlier by Radzihovsky et al. [6j [7] and 
Romans et al. jH] in the context of resonant Bose gases 
in the continuum, i.e. at low densities. This, however, 
turns out to be challenging due to the poor stability of 
the molecular Bose gas close to the resonance [9]. Here, 
we encounter a weak coupling analog of this scenario on 
the lattice, in which the stabilization of the system is pro- 
vided by the blockade mechanism leading to the 3-body 
hardcore constraint. Besides this feature, the presence 
of the lattice leads to intriguing enrichments compared 
to the continuum physics, as we will demonstrate in this 
paper. 

The qualitative picture for the Ising transition can be 



obtained within a simple Gutzwiller approach, in which 
the three-body constraint is easily built in via choice of 
the ansatz wave function [3]. However, this treatment 
leaves a number of questions unanswered, which arise on 
various length scales in the problem, and misses out im- 
portant - even qualitative - aspects of the phase diagram 
as we will show. On the microscopic scale, this concerns 
the bound state formation, as well as the correct form 
of the effective theory for dimers in the strong coupling 
limit. On the intermediate scales, relevant to the ther- 
modynamics, one may wonder to what extent the phase 
border obtained within the mean field is quantitatively 
accurate. Finally, a thorough analysis of the competition 
of the long range low energy degrees of freedom is nec- 
essary to answer the question of the true nature of the 
phase transition. We note that all these effects are tied 
to interactions, thus not available in a simple spin wave 
extension of the mean field theory. 

This paper is the second one of a sequence of two re- 
lated papers. In Ref. [11] , we have developed a quan- 
tum field theoretical framework which makes it possible 
to analytically address the above questions in two and 
three spatial dimensions. It is based on an exact map- 
ping of the constrained lattice boson model to a coupled 
theory of two unconstrained bosonic degrees with poly- 
nomial interactions. In the related paper [llj . we have 
concentrated on the formal development of this mapping, 
and performed calculations in the "vacuum limits" cor- 
responding to zero and maximum filling n = 0, 2, which 
are characterized by the absence of spontaneous symme- 
try breaking. In the present paper we apply this for- 
malism to the many-body problem. We concentrate on 
the three interaction-related aspects of the many-body 
problem mentioned above: First, we address the quan- 
titative question of shifts in the phase border, with the 
result that they are pronounced at low densities, while 
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FIG. 1. Phase diagram for the attractive 3-body hardcore 
constrained attractive Bose-Hubbard model. The black curve 
below is the mean field result, while red and blue curve corre- 
spond to d — 2, 3 where fluctuations are included. The bound 
state formation of dimers (n = 0) and di-holes (n = 2) takes 
place at the red and blue crosses for d = 2,3, determining the 
endpoints of the critical lines (Sec. Hit A bicritical point, 



characterized by energetically degenerate but different orders 
(superfluid and charge density wave) is reached asymptoti- 
cally at half filling. It can be detected experimentally ramp- 
ing a superlattice (Sec. IV I. An Ising quantum critical point, 



connecting the two ordered phases, is predicted in the vicinity 
of half filling, while the correlation length is large but finite 
away from this point (Sec. [vb. 



basically absent as the filling increases to its maximum 
n = 2. Second, making use of the perturbative results 
obtained in we consider the many-body physics in 
the strong coupling regime, and predict the existence of a 
new collective mode at half filling n = 1 , whose presence 
results from a symmetry enhancement from the conven- 
tional phase rotation symmetry U(l) ~ SO(2) exhibited 
by bosonic systems to an 5*0(3) symmetry. We propose 
an experiment to test this scenario, exploring its conse- 
quences both analytically as well as using exact numerical 
methods in one dimension. These studies lead us to call 
the system in this regime a "continuous supersolid" - a 
supersolid with a tunable ratio between the superfluid 
and the charge density wave order parameters (cf. an 
analogous phenomenon in magnetic |12[ 1 1 3 j and attrac- 
tive fermion |14j systems). Third, in a long wavelength 
analysis the phase transition turns out to be first order 
for generic fillings due to the Coleman- Weinberg mecha- 
nism [15] . This is in line with the low density continuum 
analysis, which has been carried out in detail in [TJ. In 
our constrained lattice system, however, we find that the 
radiatively induced first order transition terminates into 
a true Ising quantum critical point in the vicinity of half 
filling, which connects the two ordered phases of atomic 
and dimer superfluid. Its origin may be traced back to a 
zero crossing of the dimer compressibility together with 
a sequence of Ward identities, thus being protected by 
symmetry. An estimate of the correlation length suggests 
a broad domain of intermediate fillings 1/2 < n < 3/2 



on which the correlation length greatly exceeds the di- 
mensions of typical optical lattices, suggesting that the 
Ising quantum critical behavior could be experimentally 
observed. Our analytical approach enables us to eluci- 
date the mechanisms behind all our findings, establish- 
ing that the latter two effects are unique features of the 
three-body constraint. Our main results are summarized 
in the phase diagram presented in Fig. [T] 

The paper is organized as follows. In Sec. [Iljwe first 
review the steps that lead from the constrained theory to 
the interacting boson theory. We then prove Goldstone's 
theorem for the effective action obeying the constraint 
principle, and formulate the equation of state. In Sec. 
TIH we P ass on to the calculation of the phase border 
beyond mean field. Sec. |IV| discusses the many-body 
physics in the strong coupling limit, and in Sec. [V]we in- 
vestigate the nature of the phase transition by perform- 
ing the long wavelength limit of the effective action. Our 
conclusions are drawn in Sec. I VI I 

A summary of our results together with a closer dis- 
cussion of experimental realizations is presented in [10] . 



II. QUANTUM FIELD THEORY FOR THE 
MANY-BODY PROBLEM 



In this section we address two aspects which are partic- 
ularly relevant for the many-body physics and have not 
been discussed in pT]: The realization of Goldstone's 
theorem in our constrained model, and the equation of 
state. To prepare for this discussion and set the notation, 
we review the construction of the quantum field theory 
in Sec. |II A[ also making the paper rather self-contained. 
The reader familiar with the construction, and the reader 
who is more interested in the physics results of this work, 
may jump this section. 



A. Review of the Construction 

The starting point for our analysis is the Bose-Hubbard 
model with a three-body onsite hardcore constraint, 

H=-J^a\ aj -iiYni + |l/5jfti(n* - 1), « t3 = 0(,l) 

(id) i ' 

Here, on, a\ are the bosonic creation and annihilation op- 
erators, J is the hopping matrix element /i the chemi- 
cal potential and U the onsite interaction energy. The 
summation in the first term is performed over near- 
est neighbors. Because of the constraint, the original 
bosonic onsite Hilbert space is reduced to the three states 
\a),a = 0,1,2. 

Following Altman and Auerbach |16j . we introduce 
three operators which generate the three onsite states, 

|a)=4 l |vac) = (a!)- 1 /2 ( a t) Q | V ac), £ P) 
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from some auxiliary "vacuum" state |vac). The operators 
are not independent but obey a holonomic constraint as 
indicated above. The Hamiltonian in terms of operators 
t a reads 

H = - J [K\ W) Kf^ + 2K^Kf 1] (3) 

(id) 

i i 

where 

^(10) _ A . K {2i) _ A , - _ A , 

i — 1 i "0,11 i — 2 i 1)^5 tl a.,i — L a i^a^i- 

Note that in this representation of the constrained Hamil- 
tonian, the conventional roles of interaction and hopping 
are reversed: while the interaction enters the quadratic 
part of the Hamiltonian, the hopping term gives rise to 
effective kinematic interactions. The representation is 
therefore ideally suited in a strong coupling limit. 

In a naive Gross-Pitaevski treatment of the Hamilto- 
nian, achieved by replacing the operators with complex 
valued amplitudes t a ^ — > f a _i in Eq. ([3]), reproduces pre- 
cisely the Gutzwiller mean-field energy, i.e. a classical 
Hamiltonian field theory for spatially varying amplitudes 
f a j, where the holonomic constraint X)i /« i/«t* = 1 is 
the normalization of the onsite wave function. We now 
show how one can introduce a convenient description 
of the theory on the quantum level. To illustrate the 
method we consider the case of vanishing density n — > 
(the generalization to an arbitrary density < n < 2 
will be given below). In this limit, it is convenient to 
express the to,, operators in terms of t% t i and £ 2 ,i oper- 
ators using the constraint. Writing tg.i = | to , -i I expit^, 
we observe that the phase (pi is unphysical: it can be 
eliminated via a local redefinition of the remaining oper- 
ators, ta ti = t aii e lipi (a = 1,2). Thus, we may consider 

1/2 

to as real, and replace io.j = > Xi = 1 — — n 2 ,i 
in k\ W \ Obviously, the square roots are impracticable 
for any quantum field theory because they give rise to 
vertices of arbitrarily high order. To eliminate this prob- 

1/2 

lem, we use the fact that the matrix elements of X i 
and Xi on our subspace are the same: either 1 or 0. 
Consequently, on the subspace we may replace 

K\ W) = 4,iV^ -> tljXi, -Yj = (1 - n M - fi 2)i ),(4) 

and analogous for the hcrmitian conjugate. More for- 
mally, the replacement can be justified by noting that 
the constraint operator is a projection, Xf = X i7 and 
that the Taylor representation for a function of such an 
operator is /(X) = /(0)(1 - X) +Xf{\)^\ With this im- 
plementation of the constraint, the remaining operators 



See 1 1 1 1 for a subtlety in deriving this formula. 



ti,t2 can be treated as standard bosonic operators acting 
in a complete Hilbert space H = Yii Hi i where Hi — 
{\rii) \rrii)}, n^, m, = 0,1,... is a bosonic Hilbert space 
for "atoms" t\ and "dimers" ti at each site i : \rii) — 

(n 4 !r 1/2 (4)™*|0), and \m % ) = (m,!)' 1/2 (4) m * |0>«. The 
onsite Hilbert spaces Hi can naturally be splitted into 
a physical subspace Vi with Uj + Wj = or 1, and 
an orthogonal unphysical one Ui with rii + rrii > 1, 
%i =Vi®Ui- Important for our construction is that the 
Hamiltonian H has no matrix elements between physical 
and unphysical subspaces, (u\H\p) = (p\H \u) = 0, where 
\p) G V = YliVi and \u) € U = JJi^i, and, therefore, is 
block diagonal, H = H-p + Hjj. As a result, these sub- 
spaces do not mix during evolution, and all quantities, 
both dynamical and statistical, factorize. For example, 
for the partition function one has 

Z = Trcxp(-f3H) = Z P + Z u (5) 

= £<p| exp(-pH v )\p) + £H eM-f3H u )\u). 
M M 

Consequently, if we find a way to discriminate between 
the physical and unphysical contributions, we may indeed 
conceive the operators ti t 2 as conventional bosonic ones. 

Such a setting is provided by using the effective action 
to encode the physical information of the theory, see e.g. 
|17j . It is defined as the Legendre transform of the free 
energy W[j] — log Z[j] (we introduce a source term j — 

(jijhhjt) and use £ = (t\,ti,tl,t 2 )): 

IT ir, • (6) 

where the new variable £ = (£) is the field expectation 
value or the "classical" field. The effective action has the 
following representation in terms of a functional integral, 

exp-r[£] ^JvS^exp-S^ + 5£]+Jj T S£, j = ™(7) 

where <5£ = £ — £, and the Euclidean action S — 
It Si t\d T ti + t\d T t2 + H[ti,t 2 ]. The Hamiltonian now 
is to be interpreted as a function for classical though 
fluctuating, time dependent fields. The last identity in 
Eq. Q is the full quantum equation of motion, and the 
equilibrium situation we are interested in is specified by 
j = 0, where no mixing between the physical and the 
unphysical sector occurs. Usually, the most general form 
of the effective action is only restricted by the symme- 
tries of the microscopic theory. Since, as shown above, no 
couplings mapping from U O V are generated, we have 
identified a means to distinguish physical vs. unphysi- 
cal contributions by writing down the most general form 
for the effective action for the physical sector by directly 
excluding couplings which would violate this constraint. 

Now we generalize the procedure to arbitrary density. 
We first follow [15] but then apply our exact procedure. 
While we have so far replaced the to operator, which 
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generates the mean field vacuum state |fi) =Yli il vac ) 
we now consider a more general mean field vacuum, 



i a 

= n ( n r ° ex p(* a ^)^L) i yac ) = n b U 



(8) 



vac). 



For site independent amplitude moduli r a , these states 
allow for the description of homogeneous ground states 
with spontaneous phase symmetry breaking: If, e.g., all 
r a =/= 0, the requirement of a fixed spontaneously chosen 
overall condensate phase cj> requires the phase relation 
9 a = a.(j); this fixed phase relation is the manifestation of 
spontaneous symmetry breaking in the Fock space. We 
can now introduce a new set of operators tf a (a — 0, 1, 2) 
in which 6q creates the mean field vacuum and will be 
eliminated. Such a transformation is performed via a 
two-parameter unitary rotation, whose rotation angles 
are chosen such that the new operators fluctuate around 
the new vacuum state and do not feature expectation 
values, 



«>t _ 



with the explicit form of the rotation matrices 



Ry 




sin0/2e 2i ^ 

1 

cos 6>/2 



-sinx/2e^ 
cosx/2 



(9) 



(10) 



A finite 8(x) corresponds to a finite amplitude in |2)(|1)) 
and we will see below how these quantities are fixed via 
the Goldstone theorem. The precise relation is 

r o =cos0/2, T\ — sin (9/2 sin x/2, r 2 = sin 9/2 cos x/2. 

(11) 

At this point we can repeat the steps described above for 
the case n = in complete analogy. The constraint is 
implemented via the replacement 

&o,< ->Xt = l- bljhj - 4,^2,i, &oA-i -> Xi-(12) 

The second line is simply a rearrangement of the holo- 
nomic constraint. The resulting bosonic Hamiltonian, 
which is then quantized by means of a functional inte- 
gral, is rather complex, and we will analyze it below. 
However, it exhibits a simple structure, 



H — Egw + H$w + Hi, 



(13) 



-Egw is the Gutzwiller mean field energy and -ffgyv de- 
scribes the quadratic spin wave theory]^] The corrections 



2 A linear contribution, as naively expected in the expansion about 
the condensate, does not occur due to Goldstone's theorem, see 
below. 



to the mean field phase diagram, as well as nontrivial ef- 
fects in the deep infrared physics which we are interested 
in here, are not captured at this quadratic level. They 
are all encoded in the interaction part H lnii . 

Choosing the qualitative form of the ground state pre- 
requisites a certain knowledge about the physics of the 
system. Equipped with the right qualitative ground 
state, we can then perform quantitative calculations be- 
yond the mean field level based on our mapping. Indeed, 
Eq. ( 13 1 suggests an interpretation of our construction 



as an exact requantization procedure of the Gutzwiller 
mean field theory. This is in complete analogy to the con- 
ventional treatment of e.g. bosonic continuum systems 
with broken symmetries, where in a first step a certain 
order parameter is chosen and the theory is expanded 
around it. However, on the lattice the right choice of 
the qualitative features of the ground state might be less 
obvious. For example, spatial modulations of the order 
parameter are possible, such as exhibited by charge den- 
sity waves. This is easily incorporated in the formalism, 
and such a situation will be indeed encountered in Sec. 

m 

In sum, we have obtained the following simple result: 
supplying the most general form of the effective action 
with a constraint principle, the evaluation can proceed as 
in a standard polynomial boson theory. Similar to sym- 
metries, the restrictions on the full theory leverage over 
from the microscopic theory. Unlike symmetries, the rel- 
evance of the constraint depends on scale, being restric- 
tive on short distances, while on long distances power 
counting arguments lead to an effectively unconstrained 
though interacting spin wave theory with two degrees of 
freedom (see below). In practical computations, we can 
evaluate a theory of standard coupled bosonic fields. This 
opens up the powerful toolbox of modern quantum field 
theoretical methods for calculations in onsite constrained 
models. 



B. Goldstone's Theorem and the Constraint 
Principle 

We will derive Goldstone's theorem from a compari- 
son of the full quantum equation of motion and the effec- 
tive potential. The latter is defined as the homogeneous 
part of the effective action, obtained by inserting tem- 
porally and spatially homogeneous field configurations 
U = r[b\ 7 b 1 ,bl,b 2 }/V d+1 (V d+1 = M d /T is the quantiza- 
tion volume, M the number of lattice sites in each lattice 
direction). The possible dependences of the effective ac- 
tion and potential on the fields are strongly restricted 
by both symmetry and constraint principle; the effective 
potential is further limited by the requirement of homo- 
geneity. We will show here that the constraint leads to an 
additional U(l) invariant on which the effective potential 
may depend with no counterpart in unconstrained theo- 
ries, but we will also demonstrate that it does not break 
the validity of Goldstone's theorem - in line with the in- 
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tuition that the microscopic constraint would not affect 
the long wavelength physics too strongly. Note, however, 
that the constraint has an impact on the long wavelength 
physics, as it is indirectly responsible for the presence 
of the Ising quantum critical point close to unit filling 
n = 1 (cf. Sect. VB). Therefore, a thorough discussion 



of Goldstone's theorem seem adequate. 

Let us construct the most general dependence of the 
effective potential on the variables 61,62- For simplic- 
ity of the presentation, we focus on a spontaneously 
broken symmetry for the dimers (9 ^ 0, tt), while the 
atoms are in the normal phase (\ — 0, tt). The latter 
field can therefore be excluded from the following con- 
siderations. There are two possible terms associated to 
the original t 2 degree of freedom that might appear in 
the effective action: Either it appears as a local com- 
bination fi2 t i — £2^2,1, or as a bilocal (in general, n- 
local) combination, such that the constraint has to be 
taken into account via proper combination with e.g. 
4 i*0v While the local combination respects the £7(1) 
symmetry, the second term must appear with a con- 
jugate partner as t% jio.i^o ^ n or der to implement 
the finite density, we now apply the rotation prescription 
t 2 ,i = sbo,i + c&2,i,io,i = c&o,i — s*&2,i and subsequently 
impose the constraint 6q,i ~^ ^Q- (Here and in the follow- 



ing, we abbreviate s — sin#/2e 



cos 0/2.) Now 



we specialize to the homogeneous part of the effective ac- 
tion: We Fourier transform the operators and restrict to 
the zero frequency and momentum part of the combina- 
tions. We then find that the effective potential can be 
written as a function of two invariants, 



U(p,tf\), (14) 
p = {sX + cb\){.sX + cb 2 ) 

= s 2 + cs(Xb 2 + b\x) + (c 2 - s 2 )4& 2 - s 2 b\b u 
A = (sX + cb\){cX - sb 2 ) 

— cs + c~b' 2 X — s Xb 2 — 2csb' 2 b 2 — csb[bi, 



where X, b 2 denote the zero momentum and frequency 
components of these field expressions, and without loss of 
generality we have chosen s real. Note that neglecting the 
constraint by setting X — > 1, and considering low densi- 
ties, s — > 0/2, c — > 1, we recover the standard quadratic 
form for the condensate density from the local combina- 
tion, p = b\b 2 , with 62 = s + 62- However, the constraint 
principle requires a more complicated form of p, as well as 
the account for a second invariant A^A. In the following, 
we will be concerned with first and second derivatives of 
the effective potential with respect to 62,62: which are 
evaluated at the physical point b 2 — b\ = b\ = b\ = 0. 
Thus, we may set X — > l,b 1 ,b\ — > from the outset. 
Now we will show that the most general dependence of 
U can be further restricted. For b 2 ^{r\ b\ j(r), we intro- 



duce the basis of hermitian fields er^r), 7Ti(r), 

M T ) = 4= (o-iW + iTiW), (15) 

= (a(q) + m(-qf) , 

which as the original fields do not carry expectation val- 
ues. Here we have used the Fourier conventions 



b 2>i (r) = / e^b 2 (q), b\.{r) = / e^b^q), (16) 

J q J q 

Xi = (T,Xi), g=(w,q), ^ = y^X)- 

We calculate the local combination in terms of these op- 
erators, 



p = s 2 + (c 2 ~ s 2 )i((T 2 + 7T 2 ) + V2csa, 



(17) 



(cr = <r(q — 0),7r = 7r(<7 = 0)) and we observe that A^A, 
to the relevant quadratic order, can be written as 



A+A = S 4 + ( c 2 - s 2 )p-2( cs )V 



(18) 



Thus, the most general dependence of the effective po- 
tential on the homogeneous fields a, tt is given by 



U(p,a 2 ). 



(19) 



Now we study the mass matrix, which can be calcu- 
lated from the effective potential as the second deriva- 
tive with respect to cr, tt. In particular, the form of the 
effective potential implies for the tt mass or gap 



d 2 U{p,a 2 ) = fd 2 P u , + (dp\ 2 ull \ 

dirdn a=-n=o \3tt 2 KdnJ J a=n=cj 

(primes denote derivatives w.r.t. the invariant p) with 



(20) 



dp 
On 



0. 



d 2 p_ 
3tt 2 



s 2 . 



(21) 



To complete the derivation of Goldstone's theorem, we 
calculate the equation of motion for a from the effective 
action, but immediately specialize to the case of homo- 
geneous fields, 



ST 



horn 



dU{p,a 2 



8a 2 



(22) 



By construction we have a — as the solution of the 
equation of motion, and furthermore in the presence of 
spontaneous symmetry breaking cs =/= 0. Thus 



W = 0. 



(23) 



This simple relation indicates the presence of the gap- 
less Goldstone mode: The tt gap calculated in Eq. (20) 



vanishes due to Eq. (23 1, and since dp/dir — 0, cf. Eq. 



(21). This property is protected by the U(l) symmetry 
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of the problem. Though the form of the effective poten- 
tial is more complicated than in the continuum at low 
densities, where the effective potential depends only on 
the low density limit of the invariant p, we can explicitly 
prove Goldstonc's theorem. 



Note, that the equation of motion (22) for a also ex- 
cludes any homogeneous linear term in this field. The 
same is true for the 7r mode. Such terms would not be 



compatible with the equilibrium condition ( 22 ) . This ex- 



cludes nonzero couplings from the homogeneous terms 
~ a, 7r or 62,^2 m the effective action. Furthermore, via 



Eq. ( f 4 ) the linear terms b 2 , b 2 are connected by the con- 



straint principle to cubic terms: only the combinations 
~ b 2 X, Xb\ occur in the effective potential. Thus, com- 
bining Goldstone's theorem and the constraint principle, 
we see that the cofficients of the terms b 2 X, Xb\ must 
vanish, i.e.: 



S 3 T 



8T 



8b 2i 5b 2>i 5b 2<i hom SbliSh,- 



ST 



SbliShjSblj hom 8b[ }i 5b lA 5b 2: 



ST 



hom 



hom 



(24) 



In the presence of an atomic condensate x 7^ 0, anal- 
ogous equilibrium conditions can be derived for the b\ 
field. 

In the symmetric phases n — 0, 2(9 — 0, 7r), no distinc- 
tion between the phase and the amplitude mode appears 
and the mass matrix is degenerate. In this case, Gold- 
stone's theorem reduces to the condition for the existence 
of a dimer/di-hole bound state. 

C. Equation of State 

The equation of state is obtained as the average over 
the particle number operator N = ^2 i fii,fii — h\i + 
2n 2 ,i • Thus, after rotation we have for the particle den- 
sity 

n = (N) I M d = 2| S | 2 + (c 2 - \smb\h) + 2(b\b 2 )\ 

+ 2cls*(Xb 2 )+s{bix)} = -™. (25) 

The second equality results from the path integral rep- 
resentation of the effective action, and is due to the cou- 
pling —pN in the microscopic action. The connected two- 
point functions are given by the traces of the full Green's 
functions for b\ and b 2 . A convenient shorthand to re- 
late the connected Green's functions to its one-particle 
irreducible counterpart is £ = (b\,bi,b 2 ,b 2 ), 

(b{h) = (66) = Tr G 12 , (b\b 2 ) = = Tr G 34 , 

8 2 T 



r _ (r (2) -1\ r (2) 
^ab — )ab, 1 ab ~ XC XC ' 



(26) 



where we suppress spatial or momentum indices. Tr runs 
over these as well as over the internal (field space) indices. 



More explicit formulae will be discussed in the next sec- 
tion. Furthermore, the three-point correlation is related 
to the one-particle irreducible three-point vertex via Eq. 
(125} (cf. e.g. 



(Xb 2 



El) 

= -(b{h b 2 ) 



(b\b 2 b 2 



(27) 



E 

a.b.c 



Tr[G la G 2 h + G3 a G4h]G4 c r 



(3) 
abc 7 



.(3) 
abc 



S 3 T 



S£aS(,bSt;c 



At this point, we stress that the parameter |s| 2 in the 



equation of state ( 25 ) must not be interpreted as the con- 



densate fraction, though the formal appearance naively 
suggests such an interpretation. Instead, 2|s| 2 should be 
seen as the classical or mean field contribution to the to- 
tal particle density, and the rest of the equation is due to 
fluctuations on top of this mean field state. A standard 
interpretation of the above equation is only possible in 
the low density limits n — > 0, 2. Omitting the three-point 
correlations, Eq. ( 25 ) reduces to leading order to the fa- 



miliar form from thermodynamics in the continuum for 
9 — > 0, while taking a similar structure for 9 — > n, 

8^Q:n = 2(S9/2) 2 + (b\h) + 2(4& 2 ), (28) 

9 ps Tr : 7i = 2 - [2(<S6>/2) 2 + (b{bi) + 2(b\b 2 )]. 

In these cases, 80/2 may be interpreted as the conden- 
sate order parameter. We furthermore observe from Eq. 
(25) that around 9 = tt/2 there is a point where fluc- 



tuations are strongly suppressed compared to the mean 
field contribution due to a cancellation. For a proper 
definition of the condensate fraction in the system, we 
can use the Gutzwiller expression for the original boson 
operator expectation, (pi) = sc, however with the value 
of 9 determined from the implicit condition Eq. ( |23[ ). 
More generally, we emphasize that Eqs. ( 23|25 1 provide 
the two exact, but implicit conditions that determine the 
two parameters 9, /x A further nonzero expectation value 
for the the single atom degree of freedom, described by 
X 7^ 0, adds a further such condition analogous to Eq. 
p3l. 



Our effective action formalism is capable to describe 
the system at any finite temperature. Calculations in 
this regime are beyond the scope of this paper, but let us 
sketch how the high temperature disordered phase is de- 
scribed within our theory. Increasing the temperature in 
the system will populate the connected parts of Eq. ( 25 ) 



increasingly such that at the phase transition to the sym- 
metric phase without symmetry breaking fx, f 2 — > 0. In 
other words, the condensate angles vanish, 6, x 0. We 
may interpret this scenario as the complete population of 
the "vacuum amplitude" /o, which is needed to fulfill the 
holonomic constraint but does not enter the equation of 
state. The effect of destruction of the order parameters 
can also be seen from the condition m 2 = d 2 U jdix 1 = 0. 
A finite temperature will act to generate a positive ther- 
mal mass or gap contribution, such that at some temper- 
ature there exists no finite 9, x an( i a gapless mode ceases 
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to exist. At this point, where Goldstones's theorem can 
no longer be satisfied, the symmetry broken phases be- 
come unstable and the system enters the disordered high 
temperature phase. 



III. ASF-DSF PHASE BORDER 

In this section we embark the calculation of the phase 
border. We will study the phase border by approach- 
ing it from the dimer superfluid side where there is no 
atomic condensate, and calculate at which interaction 
strength the atoms become unstable towards an atomic 
superfluid. Thus, we first provide the explicit form of 
the Hamiltonian in the presence of a dimer superfluid, 
but for atoms in the normal phase. We then consider the 
low density limits n — ¥ 0, 2. In these limits, we can estab- 
lish a controlled small density expansion describing the 
deviation from exactly n = 0, 2. The central objects for 
the discussion are the atomic and dimer (di-hole) Green 
functions, which we know exactly in the limits n = 0, 2 
[TT] . The analysis reveals the intuitive result that the 
leading many-body effect is a modification of the vac- 
uum (n = 0, 2) Green functions due to the condensate 
mean field. The dominant fluctuations in these limits are 
thus vacuum fluctuations renormalizing the Green func- 
tions, while the many-body effects can be captured in 
terms of a Bogoliubov or spin wave theory. More specif- 
ically, we find that vacuum fluctuations strongly modify 
the relation fx(U) compared to the mean field relation 
fj,(U) = —U/2, while the role of many-body effects con- 
sists mainly in depletion effects in the equation of state. 



We find that the high energy vacuum fluctuations have a 
much more pronounced quantitative effect on the phase 
border than the condensate depletion in the limit n — » 0. 
For n — > 2 instead, both effects are rather small, which 
may be understood in terms of an already tightly bound 
di-hole state in the region of atom criticality. Based on 
these insights, we do not expect a strong shift in the 
phase border in the region n = 1, which takes place at 
even stronger coupling, and thus more deeply bound two- 
particle states. We therefore propose an extrapolation of 
the scheme from the controlled limits n sa 0, 2 to the 
intermediate regime nwl. 



A. Rotated Hamiltonian for the dimer superfluid 
phase 

Let us now focus on the phase border to the dimer su- 
perfluid state. As anticipated above, we address it from 
the DSF side where |1) is not macroscopically populated 
and thus % = 0. The kinetic and potential energy oper- 
ators read, in the new field coordinates, 

K\ XQ )= bUcb 0ti - s*b 2<i ), Jf< 21 > = (s%. + cbUb^, 



R 



(n)_ 



b\ ,6 



E 



(22) 



( s * b 0,i + cb 2,i)( sb 0,i+ cb 2,i), 

(29) 

with c = cos#/2, s = sin#/2e 21 ^ as above. We can now 
write the Hamiltonian operator in terms of the new vari- 
ables, and implement the constraint via &o,i ~^ = 
(1 — fii i — n,2,i), absorbing the phase of bg into the re- 
maining two degrees of freedom as discussed in Sec. |II A| 
Further making use of the projective property X 2 = Xi 
we find 



H[b 1 ,b 2 ]=H^'+H ] 



(10) , rr(21) 



kin ^k 



(split) 



pot ; 



(30) 



#ki„ 0) + H^ ] =-J £[(c 2 + 2\s\ 2 )bl l X l X J b l . J + (2|.s| 2 + r 2 - 3c (sb^X^b^ + s^b^Xjh, 



{id) 



#£piit) = -^2jJ2 [(c 2 - H^/mX^ij + c (s*X i b 1>i X j b ltj - sbl ti b lti bl tj b ld ) + h.c] , 

{hi) 

tfpot - (-2M + U)M\s\ 2 + [-2n + U)(c 2 - \s\ 2 ) ]T 4/2,, + [-H - (-2/i + U)\s\ 2 )} ^ b\j x . 

i i 

+(-2/i + U)cJ2 [sbl,iXi + s*X l b 2 ^ . 



We remind the reader that, as shown in Ref. and and the cubic parts are connected to the linear ones via 



briefly discussed in Sect. |II A these operators b\ t 2 may 
be interpreted as standard bosonic operators. The cubic 
term in the second line of -ff po t can be omitted from the 
outset, and we will do this in the following: As argued 
above, the coefficient of the linear part has to vanish due 
to the equation of motion, i.e. the equilibrium condition, 



Eq. ( 24 ) , such that their coefficient has to vanish as well. 



However, this does not exclude the possibility of nonlocal 
cubic terms as they appear in the kinetic terms of the 
Hamiltonian. The total Euclidean action in the presence 
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of condensation reads 



S[h,b 2 } = / dT(^2b[ ti a r bx,i + ^,idrb 2 , i + H[b 1) ba] 



(31) 



where the Hamiltonian is to be interpreted in the Heisen- 
berg picture and as a function of classical field variables. 
Quantizing this theory with the functional integral leads 
precisely to the representation of the effective action Eq. 

0- 



B. Low density limits n « 0, 2 

In the following, we will analyze the theory in the vicin- 
ity of the physical vacua where n sa 0,2, described by 
6>o ~ 0, 7T. The limits n = 0, 2 have been discussed in [TT] 
in detail. The Hamiltonians governing these situations 



describe the scattering of few particles in the absence of 
many-body effects and can be written as 



H n = 



[9n,lK ti XiXjblj + .9„,24,Aj &t M & 2,' 

(id) 

V2J(b\ ib^Xjtu + b\ .Xjb\ 



(32) 



The operators b\^ represent the bosonic single and two- 
particle excitations, corresponding to atoms and dimers 
resp. holes and di-holes. Here, for n — we have go.i — 
J, 50,2 = 2J,^ = M and for n = 2, g 2 ,i = 2 J, g 2< 2 
./. fj, 2 = —fi + U. At these points the exact solution of 
the (two-body) scattering problem, and thereby an exact 
calculation of the atomic and dimer Green's function, is 
available as shown in [11] . While the case of the atomic 
Green function is trivial as there are no renormalization 
effects in the vacuum, for the dimers/diholes we find the 
results 



G~ l 1 (u>; f M>,k) = U + 



d d q 



(27r) d -2(e q 
G^V,M2,k) = U+ ^(kj-2/i 2 ) + i 



d d q 
(27r) d 



2flQ 



-4(e q + e q _ k ) 



2/x 2 



(33) 



We will now perform a controlled expansion in the con- 
densate angle deviation 59 <C 1 from the special points 
#o = 0, tt. It corresponds to a Bogoliubov approximation 
for the condensation physics, but with coupling constants 
obtained from the exact solution of the "vacuum" scatter- 
ing problem. The procedure amounts to a resummation 
of ladder diagrams. These vacuum fluctuations are re- 
sponsible for strong shifts in the phase border as we will 
see. 

Our expansion is defined with Hamiltonians of the form 



H — H„ 



2 SH « 



o{se 2 ). 



(34) 



The additional Hamiltonian 5H n generates new scatter- 
ing vertices which are 0{89). Diagrams with more than 
one of the new vertices may thus be discarded, and we 
may restrict our attention to diagrams with at most one 
of them. They will be discussed in a moment. 

In the low density cases n « 0, 2 the discussion can be 
lead in parallel, due to the similar mathematical struc- 
ture of the Hamiltonians. Around #o = 0, we replace 
c = l,a = 59/2 and the additional Hamiltonian reads 

SH Q = -J^V^lXib^Xjhj - blb.Xb^: 

(id) 

+X i b 1<i bl j b 2d + h.c. (35) 



r 



Similarly, around 9q = tt, setting s 
find 

SH 2 = ~SH . 



= l,c = -89/2 we 



(36) 



Note, that the zero order Hamiltonians H n are related 
by a more complicated transformation of parameters, re- 
flecting the absence of a particle-hole symmetry. 

Let us now discuss the impact of the additional Hamil- 
tonian. The stability of the ASF phase is encoded in the 
full atomic mass matrix, i.e. the inverse Green's func- 
tion Gi (u>; fx, q) at zero frequency and momentum: If 
the eigenvalues of the mass matrix are all positive, the 
phase without atomic condensate (i.e. the condensed 
dimer phase) is stable. The instability towards a state 
with atomic condensate can thus been inferred from the 
vanishing of an eigenvalue of this matrix, or 



detG^V = 0;^,q = 0) = 0. 



(37) 



Thus we discuss the beyond mean field effects modifying 
the inverse atomic Green's function. From the exact so- 
lutions of the vacuum problems at #o = 0, tt we know that 
in these limits the inverse atom Green's function is not di- 
rectly renormalized: there are no diagrams in the vacuum 
limits which cause renormalization, but clearly, the func- 
tion /i(U) entering the atom propagator changes when 
taking the exact dimer or di-hole Green's function into 
account. We now concentrate on the effects of 8H n . At 
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linear order in 50, we find a direct condensate contribu- 
tion to the atom inverse propagator on the off-diagonal. 
This is the contribution familiar from Bogoliubov theory 
in the low density limit, and we see that our generaliza- 
tion to arbitrary density produces such a structure also 
at high density. Now we have to consider the effect of the 
new vertices. As argued above, we can restrict ourselves 
to diagrams carrying a single one of them. We focus on 
diagrams which renormalize the inverse atom propaga- 
tor. These are tadpole diagrams. The diagrams renor- 
malizing the diagonal entries must be 0(S9 2 ) in order 



to ensure particle number conservation. The diagrams 
renormalizing the off-diagonal entries must involve one 
of the new vertices 0(50), and the trace over the inner 
line scales with a function f(89) with /(0) = 0. Con- 
sequently, the fluctuation contributions are higher than 
linear order for both diagonal and off-diagonal entries 
and can be discarded. Thus, the full atomic mass matrix 
at O(50) reads (we separate true potential (binding) en- 
ergy from kinetic energy, ^(U) = Eb(U)/2 — Jz at n « 0, 
fi(U) = -E b (U)/2 + U + 2Jz at n w 2 [IT], z = 2d the 
lattice coordination number) 



C- 1 (u> - 0- u a - m - ( -^( f/ )/ 2 ^V2Jz50/2 \ 
r- 1 (u o- u a m - ( - E ^ u )l 2 -1V2JZ50/2 \ 

G 1A =>-°^' ( l- )-( v _2^/2 -E b (U)/2 ) 



(38) 



Hence we conclude that the dominant effect beyond mean 
field theory, which implies the simple linear relation 
n(U) — U/2 for the binding energy (see but the argu- 
ment is repeated below Eq. ( |43[ ) for convenience), comes 
from vacuum fluctuations, which determine the value of 
n(U) (or, equivalently, E b (U)) as a function of the in- 
teraction strength U. These high energy fluctuations are 
responsible for shifts of the critical point. Note that the 
sign on the off-diagonal of the second equation is unphys- 
ical as it can be absorbed by a phase rotation of the order 
parameter. Thus, the equations have the same form. 

The critical interaction strength can now be extracted 
from the characteristic equation ( 37 ) which reads explic- 
itly, in both limits, 



(~E b (U c )/2) 2 - 2(Jz6 



0. 



(39) 



The physical solution is given by E b (U c ) = — 2^/2 Jz50, 
i.e. in the vacuum limit 50 — > also the binding energy 
vanishes. Thus, the critical interaction strength for the 
formation of the dimer or di-hole bound state coincides 
with the energy scale of the single particle excitations 
(atoms or holes) becoming critical. In these limits, we 
may quantitatively estimate the dependence of the in- 
teraction strength on e.g. the condensate fraction in two 
and three dimensions. In d — 3, the binding energy starts 
quadratically due to the non-analyticity in the fluctua- 
tion integral. In contrast, in d = 2, the fluctuation inte- 
gral features the well known logarithmic behavior. This 
yields, in dimensionless units, the physical solutions 



n^0: d = 3:U c = Uo(l + 2 3/i crv r 5d), 

d = 2:U c = 2n(\og^y 1 

n^2: d = 3:U c = U 2 (l + 2- 3 / 4 cr%=^\ 
v |f 2 | 

rf = 2:[/ c = -3 + 7r(log^)' 1 



(40) 



U/Jz 



Atomic Superfluid 





FIG. 2. ASF-DSF phase boundary for the attractive 3-body 
hardcore constrained Bose-Hubbard model. The black lower 
curve is the mean field result, while red upper and blue mid- 
dle curve correspond to d = 2, 3. The crosses denote the 
endpoints of the critical lines. For high densities, the d = 2, 3 
results are very close to each other except for a tiny region 
close to maximal filling. 



with numbers Uq w -4/3, a w 0.42, A w 5.50, and U 2 ~ 
— 11/3 [TT]. We note that due to the formation of the di- 
hole bound state at a finite interaction strength despite 
the logarithmic dependence of the fluctuation integral, 
also the critical interaction strength for n — > 2 remains 
finite. 

Importantly, we find a non-analytic dependence of the 
critical interaction strength on the condensate density 
50 which is due to the strong fluctuation effects in the 
vicinity of the bound state formation. This is in contrast 
to the mean field result, which shows a linear depend ence 
on the condensate angle 50 = ^/2n (see [3] and Eq. ( 45 ) 
below). 
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C. Phase Diagram 

The analysis of the low density limits in the last para- 
graph reveals that the strongest beyond mean field effect 
comes from the high energy vacuum fluctuations renor- 
malizing the relation (i(U) away from its mean field value 
fx = U /2. In our practical implementation of the calcula- 
tion of the phase boundary, we will rely on this separation 
of vacuum and many-body effects. For the calculation of 
the phase diagram at a fixed total density n, we now also 
take modifications of the mean field equation of state 
n = 2 1 s | 2 into account. We discuss the equation of state 



(251 in an approximation where we omit the three-point 



correlations from the outset. Furthermore, in our con- 
crete computations, we restrict ourselves to the calcula- 
tion of the atomic depletion (b\bi). We will find that this 
contribution is small compared to the condensate part 
and thus has a small influence on the phase border only 
- the dominant effect shifting the border stems from the 
vacuum fluctuations leading to a strong modification of 
fi(U). Based on this observation, we do not expect that 
the dimer depletion strongly modifies the phase border. 

In order to calculate the atomic depletion, we first con- 
sider the quadratic spin wave theory, 



Si[b\M 



(hi) 1 



(b\( q )M-q)) 



A, 

i q -iw + p q 

„2\ Ttjl 



b\(-q) 



Pq=Po+5p q , Po = - f i(l-2s z )-Us z - Jz(l + s*), 5 P(l = 2(1 + s 2 )<Se q , 
A q = A + (5A q , A = 2y/2Jzcs, <5A q = 4v / 2cs£e q , 
(5e q = J ^(1 - cosg A )- 



The matrix in the second line is the inverse atomic Green 
function in frequency and momentum space G 1 " 1 (w; //, q). 
With these preparations, the approximate equation of 
state and the atomic depletion is found to be 



n = 2| S | 2 + (c-| S | 2 )(6^ 1 ), 



P q -A q 



(42) 



where we have performed the frequency integral by clos- 
ing the contour in the upper half plane. 

As stated above, we consider the renormalization ef- 
fects on the inverse atom propagator which are present 
already in the vacuum problem. These are encoded in 
the value of the chemical potential fJ.(U) and depend on 
dimension. The condition determining [i reads 



G- / 1 h (cj = 0;M ! q=0)=0, 



(43) 



where in the vicinity of n = we use the full dimer 
Green's function Gd, and close to n = 2 the di-hole ex- 
pression Gh as given in Eq. (33 1. In contrast, in the 
mean field approximation which uses the "bare" inverse 
dimer propagator iu — 2/i + U, the above condition eval- 
uates to n{U) = —\U\/2 independent of dimension and 
of whether we are close to zero or maximum filling. The 
critical point is determined by the atoms becoming un- 
stable towards condensation. This is indicated by the 
condition 

detGr 1 (^ = 0; M ,q = 0)=p 2 -A 2 = 0. (44) 



We solve the system of equations ( 42|43|44 1 in two and 
three dimensions numerically. In particular, Eq. (43) de- 



couples from the other two equations within our approx- 
imation scheme yielding the renormalizcd relation n(U), 
such that we merely need to solve ( |42|44 ) with fx(U) as 
an input. The result for the phase diagram is plotted in 
Fig. [TJ We compare these results to the mean field ap- 
proximation, which uses \i = U/2 and n = 2|s| 2 for the 
equation of state, resulting in the critical interaction 



Uc 

Jz 



= -2(yJl-n/2 + Jnf 



(45) 



As anticipated above, the beyond mean field effects are 
mainly due to the fluctuations accompanying the bound 
state formation, which strongly modify the relation fi(U) 
as compared to the mean field value, while we find a sub- 
dominant role of the many-body depletion effects. The 
shape of the phase boundary directly reflects the non- 
analytic, dimension dependent behavior associated to the 
bound state formation. The quantitative effect is more 
pronounced below half filling than above. This may be 
traced back to the fact that the domain where nonpertur- 
bative fluctuation effects play a role is smaller in the high 
density regime than at small densities, cf. [11 . A simple 
picture can be given as follows: In the limits n — > 0, 2, 
the criteria for the atom criticality (zero eigenvalue of 
G^ 1 ) and the microscopic bound state formation (zero 
eigenvalue of G^ h ) coincide and fluctuation effects on 
the phase boundary are substantial. Moving away from 
these limits, the absolute value of the critical interac- 
tion increases, and the microscopic bound state is al- 
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-0.04 - 



0.0 0.5 1.0 1.5 2.0 

FIG. 3. Condensate depletion due to the atomic two-point 
function in two dimensions. The overall shape of the curve 
with zero crossing and sign change is determined by the func- 
tion c 2 — |s| 2 premultiplying the two-point function, cf. Eq. 
( ]25[ l. The small overall size shows that the depletion effects 
produce only very tiny corrections to the phase border. 



ready tightly bound at the point where atom criticality 
is reached. The critical line then approaches the mean 
field phase boundary up to small perturbative corrections 
Thus, though our approximation is lacking a strict order- 
ing principle when moving away from the limits n — > 0, 2, 
we therefore expect the mean field result to be rather ac- 
curate. 



IV. MANY-BODY PHYSICS IN THE STRONG 
COUPLING REGIME AND A CONTINUOUS 
SUPERSOLID 

In this section we investigate the system in the strongly 
correlated limit J/\U\ — > 0. In particular, we identify a 
bicritical point at half filling of atoms (n = 1), at which 
homogeneous superfluid order (spontaneous phase sym- 
metry breaking) and charge density wave order (spon- 
taneous translation symmetry breaking) are degenerate. 
The bicritical point is due to a symmetry enhancement 
from the conventional U(l) ~ SO(2) to SO(3), which 
is seen to be intimately connected to the 3-body con- 
straint. We term the system in this regime a "contin- 
uous supersolid" due to the degeneracy of phase and 
translation symmetry breaking orders, where the order 
parameter may be rotated continuously from one to the 
other without energetic cost. This behavior is in contrast 
to other occurrences of supersolidity in bosonic systems 
[ITJ] . Though this state is only reached asymptotically, it 
governs the physics in strong coupling and close to half 
filling, and we work out the observable consequences of 
this situation. We also propose a simple experiment to 
verify this scenario. 



A. Analytical Approach 

Before embarking the calculation, let us stress that 
our beyond mean field approach is indispensable to set- 
tle these issues. Indeed, a straightforward compari- 
son of the simple Gutzwiller mean field energies of the 
dimer superfluid and a charge density wave state (CDW) 
yields degenerate energies for these two states for all fill- 
ings: The Gutzwiller mean field CDW state is given by 
| CDW) = Iljcvcn l 2 )i|°)i+i) and for a nxed average den- 
sity n has energy density BcDw/^ d = (U/2)n, which 
precisely equals the mean field energy density Eusf/M 4 
of the dimer superfluid for all particle densities. In con- 
sequence, the question of the correct ground state cannot 
be decided within the simple Gutzwiller mean field theory 
(though a superfluid is clearly more plausible for incom- 
mensurate fillings). It is necessary to first integrate out 
the high energy single particle degrees of freedom, mak- 
ing the dimers true propagating and interacting physical 
excitations. Moreover, even the second order perturba- 
tion theory is not fully conclusive as we will see. The 
deviation from the second order result, calculated in [llj . 
plays a key role in the following discussion. 

In [11] we have calculated the effective theory in the 
strong coupling limit: Perturbatively integrating out the 
single particle excitations up to fourth order, and tak- 
ing the constraint principle for the effective action into 
account, we found the low energy effective Hamiltonian 

H ctf = -t ^ (tl ti XiXjt2,j - \fi 2 ,ifi2,j) - V-d ^ n 2 ,i 
(id) 

(46) 

with t = 2J 2 /\U\, fid the effective dimer chemical po- 
tential and the dimensionless ratio of nearest-neighbour 
interaction to hopping A = v/(2t) discussed below. Since 
in the perturbative limit there arc only virtual single par- 
ticle excitations, we may replace the constraint operator 
Xi = 1 — nij — n2,j — > 1 — «2,i- In this case we have the 
following mapping to effective spin 1/2 degrees of free- 
dom, which will more clearly reveal the physics of the 
model, 

sf = i**+ify = (-ytt ij X j , (47) 

s z j = n 2tj - 1/2, 

where on the bipartite lattice with sublattices A and B 
we use (— y = + for j 6 A and (— ) J = — for j 6 B. Up 
to a constant the Hamiltonian then takes the form 

H cS =2tJ2 + s\sV + Asfs*) fi d 4- (48) 

(iJ) 1 

The anisotropy parameter A = v / (2t) evaluates to A = 1 
in the second order perturbation theory, corresponding to 
an isotropic antiferromagnetic Heisenberg model - note 
the sign change in t due to the sublattice dependent sign 
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in ( 47 ) . The fourth order calculation yields [^] 



A = l-8(z-l)(^) 2 <l. 



(49) 



This result has been derived in [TT], cf. Sect. V.D.2, 
Eq. (60). It is obtained as the ratio of dimer-dimer 
interaction and dimer hopping coefficient calculated at 
fourth order. Next-to-nearest neighbour terms also ap- 
pear at fourth order, but their numerical coefficient is 
much smaller than for the nearest neighbours due to the 
restricted pathways contributing to these terms, and are 
thus neglected. 

For A = 1 and half filling of atoms n — 1, where the 
term involving the chemical potential vanishes, the sys- 
tem exhibits a symmetry enhancement from 5*0(2) ~ 
U{\) (corresponding to rotations in the x — y— plane, or 
phase rotations, generated by exp i9 z S z cx exp i9 z N, with 
global operators S a = Ei s ?>^ = E>2,*) to SO(3) 
(corresponding to arbitrary rotations on the Bloch sphere 
expiJ2 a 9 a S a ). 



The S f O(3)-invariance of the quadratic part in (48) 
also implies a simple transformation behavior of the to- 
tal Hamiltonian under a discrete particle-hole or charge 
conjugation transformation: The special choice U = 
exp(inS x ) implements the mapping 



(50) 



Under such a transformation TV — > M d — N and, hence, 
S z —> M d /2 — S z . The particle-hole symmetry makes 
the phase diagram of deeply bound dimers symmetric 
under the replacement fi2 — > 1 — fi2- In general, such a 
symmetry is absent. Moreover, it is also not present in 
the opposite limit of strong repulsive interactions, which 
is asymmetric when fi\ is replaced with 2 — hi- 

The S'0(3)-invariance is a peculiar feature of the lead- 
ing order perturbation theory. Its physical origin is well 
understood in terms of the geometric argument which 
relates hopping and interaction paths, cf. Sec. IV D in 
the companion paper |11) . At second order, no other in- 
teraction processes can occur, and thus the hopping and 
interaction constants must be equal. However, as seen 
in [TT], at fourth order perturbation theory additional 
interaction processes yield A < 1, thus reducing the sym- 
metry to SO (2), and also spoiling the particle-hole sym- 
metry. In addition, several other terms are generated, 
which describe next-to-nearest neighbour hopping and 



3 One may wonder about implicit density effects for the perturba- 
tive calculation, that is, if the perturbative calculation at zero 
density is sufficient for the calculation of the effective theory 
for all densities. This may be discussed by studying the limits 
n = 0, 2. At second order, the emergent particle-hole symme- 
try ensures that the perturbative results coincide for both cases, 
while at fourth order, differences occur, pointing at the above 
mentioned implicit density effects. However, due to the identical 
diagrammatic structure one still finds A < 1 for n = 2, which is 
the crucial ingredient for the argument presented here. 



interaction, or three- and four-spin interactions. Never- 
theless, the proximity to the Heisenberg point A = 1 has 
an impact on the phase diagram, and we will use its well 
know properties to understand the phase diagram and 
the nature of the low lying excitations in the perturba- 
tive regime. 

The order parameter for this model is given by the ex- 
pectation value of the Neel vector JV a = 2~2j{— ) 3 sf. Its 

vector character is under 5*0(3), [S" 1 , A/'' 3 ] = lea^N 1 , 
i.e. global spin rotations transform the Neel vector com- 
ponents into each other. Translating back to the origi- 
nal boson language, a finite (N z ) corresponds to charge 
density wave order, while finite values of (M x ), {N v ) in- 
dicate dimer superfluid (DSF) order. For the isotropic 
Heisenberg model without magnetic field (or at half fill- 
ing), [H,S a ] = for all a, and thus CDW and DSF or- 
der are degenerate. The perturbative limit of our model 
thus realizes a bicritical point 20J with two competing or- 
ders. Such an enhancement of internal symmetries is well 
known in magnetic systems [121 113] , but less common and 
intuitive in physically realizable bosonic models, which 
usually only exhibit phase symmetry. Due to the degen- 
eracy of phase and translation symmetry broken states, 
both order parameters are generically nonzero, and we 
may term the state a continuous supersolid, whose ex- 
perimental implications are studied below. 

We can make this discussion even more explicit when 
changing from the spin to a hardcore boson language 

87 = (-1)^, sj = (-lM 4 - h% - i (51) 



f 2 



0. We con- 



where the hardcore bosons obey h 
sider infinitesimal SO(3) transformation with the pa- 
rameters e a <C l(a — (x,y,z)), Ss" — i [S 1 , s"] , where 

S = YliB £ P s i ■ The explicit form of the above transfor- 
mation reads (s • = s x ± is^) 



5s] 



IE S j %£ z Sj ; 

ie*s z + ie z sj , 
_i e * s + 4. i £s - 



(52) 



with e = e x + is y . In terms of bosonic hardcore operators 
one thus obtains 



8h 



Shi 



(-iy 
l-iy 



-ie(hthj 
ie*{h%- 



- 1 £*h' 



F 

2' 



4 



ie z hj 



(53) 



Note that the last terms in the first and second lines are 
just usual gauge transformations. 

The change in the operators results in the change of 
their mean-field values. If one introduces the usual super- 
fluid order parameter -0 = ^ . (hj) and the CDW order 

parameter J\f = (J\f z ) = E^- 1 )^ (( h \ h j - then the 
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corresponding change in the order parameters is 



Stjj = ie z ip + ie*Af, 
SJV = ~\e*^* + ietf. 



(54) 



It is easy to check that the above transformation leaves 
the combination M 2 + \tp\ 2 invariant. We therefore can 
conclude that the SO (3) symmetry corresponds to canon- 
ical transformations of the dimer operators, which change 
both superfluid and CDW order parameters, but leave 
the combination Af 2 + invariant. 

Another important point about the 50(3) symmetry 
is that it is broken for n ^ 1 not on the Hamiltonian level 
(the Hamiltonian with A — 1 is always 50(3) symmetric) 
but on the level of the subset of states (with a fixed S z ), 
on which it has to be minimized. In a generic case n 1 
(and, hence, S z 0) the subspace reduces the symme- 
try down to 50(2) ~ U(l) gauge group. In the case 
n = 1 with S z — 0, however, the subspace contains the 
manifold of spin-singlet (and, therefore, 50(3) symmet- 
ric) states, which have the lowest energy. The symmetry 
transformation corresponds simply to the motion on this 
manifold. 

A similar scenario (an enhancement to a pseudo SU (2) 
symmetry) is actually observed in attractive lattice 
fermion systems [14] . Similar to the fermion system, the 
symmetry enhancement is thus a unique consequence of 
the 3-body hardcore constraint. Indeed, attractive lat- 
tice bosons without such constraint, analyzed in detail 
by Petrosyan, Schmidt et al. [ST], show a different be- 
havior: Due to the possibility of virtually occupying a 
lattice site with three atoms, it is found A = 4. This 
places the unconstrained attractive bosons in the "Ising 
limit" , which was analyzed further in the latter reference. 

As we find A < 1 in fourth order perturbation theory, 
the bicritical point is approached from the homogeneous 
superfluid, which is energetically favoured over the charge 
density wave. Nevertheless, we may expect important ob- 
servable consequences. Indeed, the symmetry enhance- 
ment 50(2) — ► 50(3) implies the emergence of a second 
gapless, and therefore collective, Goldstone mode. For a 
weakly explicitly broken 50(3) symmetry, one still has a 
near gapless collective mode with experimentally observ- 



able consequences. We propose an experiment, which is 
based on the idea of explicitly rotating the macroscopic 
Neel vector from the x — y— plane representing superfluid- 
ity to the z— axis, realizing a CDW ordered state. We will 
also show that this experiment allows to quantitatively 
characterize the pseudo Goldstone mode. 

To favor CDW ordering, we explicitly break the lattice 
translation symmetry via introduction of a superlattice 
shifting the single particle energies on adjacent sites: 



n 2 ,i -> 



-fJLA ^ fl2,j 



[IB fl2,j 



(55) 



-/'- 



E 



n-2,: 



n 2 ,; 



[i-A 



jGB 



2 .2 

Here, fi is and average chemical potenial and v an imbal- 
ance parameter. Now we calculate the mean field ground 
state as well as the spectrum of excitations of the effective 



low energy Hamiltonian ( 46 1 , using the rotation formal- 
ism (cf. Sec. II A). The approach is fully equivalent to 



the leading order 1/5 expansion, which is not a well con- 
trolled expansion for 5 = 1/2, but is known to yield the 
main features of the Heisenberg model in external fields. 

At the Heisenberg point, the 50(3) spin symmetry 
requires an enlarged parameter space for the order pa- 
rameter describing the ground state of the system. We 
consider an ansatz parameterized by two angles for the 
rotation of the t 2 degree of freedom (in contrast to the 
U(l) ~ 50(2) case with a single rotation angle for t 2 ), 



cos(0 + <#)/2 sin(6l + ( < 5 ; )/2e 2i ^ 
-sin((9 + ipi)/2e- 2 ^ cos(0 + y>j)/2 



(56) 



Here, I = A, B is an index which depends on the sublat- 
tice A or B, thus enabeling the description of a spatially 
modulated phase. For example, the choice — 0, tpA — 
0, ipB — 7T describes a charge density wave. The homoge- 
neous choice 6 ^ 0, tpA = ^5=0 describes a superfluid 
ground state. The rotation matrix is only 2x2 since we 
have integrated out the atoms. 

Expanding the thus rotated Hamiltonian, and replac- 
ing bo t i — > Xi = 1 — fi2,i, to second order we obtain 
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E mi /M d = -2tz(s A c A s B CB - A44) - |(4 + s|) + -(si - 4), (57) 
Hlin = [tz( ~ c B s B {c 2 A - s 2 A ) + 2Xs 2 B s A c A ) + (-fi + D)c A s A ] ^(4,i + b 2,i) 

+ [tz( - c A s A (c 2 B - 4) + 2\s A s B c B ) + (-(i ~ v)c B s B ] y^X b l,i + b 2j), 
H sw = o X] [(~ ( c a c b + s a s e) - 1\c A s A c B SB){bl tj b 2 ^ + 4,i & 2,j) 

<i.i> 

+ ( C A S I + C I S A _ 2Xc A S A C B S B )(b 2 ,jb 2 ,i + &2,i4,j)] 

+ (2iz(2cASyiCBSs - As|(4 - 4)) + (-/■* + ^)(4 - 4)) X! & 2,i & 2,i 
+ (2iz(2cASyiCBSs - As^(4 - 4)) + ( — A* ~ D )( C B ~ s %)) X! & L & 2,i- 



Here s; = sm(9 + (pi) etc., and we have set the sponta- 
neously chosen phase </> = without loss of generality. 
M is the total number of sites in each lattice direction. 
Since we break the lattice translation symmetry via our 
choice of the ansatz for the ground state, it is important 
to be careful with the position indices - 62, i are located 
on the sublattice A, b^.j on the sublattice B. Eventually 
we are interested in a situation at hxed density. The local 
density operator to quadratic order takes the form 

4,i*a,< = *» + ClSl ( b h + 62 .' ; ) + ( c ? ~ s ?)™2,i, (58) 

where I = A(B) for i £ A(B). Thus, in the mean field 
approximation the equation of state reads 

n = N/M d = |(4 + 4)- (59) 

At half filling n = 1, this implies s A = c B ,s B = c A . 
Together with the relations sf + cf — 1, within this ap- 
proximation everything may be expressed in terms of e.g. 
s A alone. In particular, the mean field energy determin- 
ing the ground state takes the form 

E mi /M d = -2tzs\{\ - 4)(1 - A) + v{2s\ - 1).(60) 

The ground state is found from identifying the stable 
minima with respect to variation in 6, and thus we have 



I 

dropped the contribution from the chemical potential, 
as it contributes a rotation angle independent constant 
only for effectively fixed density. For v — 0, the ground 
state for A < 1 is the homogeneous superfluid with s\ = 
1/2. For A > 1, the charge density wave with s A = 
0, s B = 1 is favored. At the Heisenberg point A = 1, both 
states are degenerate in accord with the exact symmetry 
argument. Now we consider the relevant case A < 1. 
Tuning v away from zero by ramping the superlattice, 
the superfluid acquires a spatial modulation, s\ = (1 + 
i/)/2, s B = (1 — v)/2 1 v = v/tz{X — 1), where at a critical 
value 

\v c \ = 1, \v e \ = tz(l - A) « 8tz(z - l)( J/U) 2 (61) 
the SF is destroyed in favor of the CDW. Thus, ramping 
the superlattice corresponds to rotating the Neel order 
parameter. As we will see below, the critical value cor- 
responds precisely to the gap of the pseudo-Goldstone 
mode. Hence, via measurement of the SF correlations 
[22] . which cease to exist at v Cl one can quantitatively 
determine the characteristic property of the additional 
collective mode. 

The chemical potential fj, is determined from the equi- 
librium condition that the linear terms vanish, evaluat- 
ing to fi = Xtz independent of v. Inserting this and the 
above expression for s A , and switching to the Lagrange 
formalism, we obtain the Gaussian action 



<7 q = tz(^(l-^ 2 ) + l)e q , /i q (u;) =kj + iz(l + ^(l-j/ 2 )? q ), e q =±^cosq A . 



The spectrum of excitations may be computed from the vanish. We obtain 
condition that the determinant of the fluctuation matrix 

w(q) = ±tz(((l + [A — 1][1 — is 2 ])^ + 1)(1 - e q )) V (b) 
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For v = and A < 1, the dispersion simplifies to 

l li 

w(q) = ±tz((\£q + 1)(1 — e q )) , and there is a single 
Goldstone mode at q = 0, corresponding to the spon- 
taneously broken E/(l) ~ 5*0(2) symmetry in the dimer 
superfluid. For A < 1, there is a second near gapless mode 
with gap tz(l — A) located at the edges of the Brillouin 
zone q = 7r. At the Heisenberg point A = 1, this gap 
closes, and the system features the two Goldstone modes 
corresponding to the spontaneously broken 5*0(3) sym- 
metry. The system is then at the bicritical point where 
the order parameter can be freely rotated on the Bloch 
sphere. In the general case, the gap of the second near 
gapless mode is given by 



TABLE I. Low energy dispersions in the perturbative regime 
(second order), z is the dynamic exponent, uj ~ |q| z . 



A = tz(l - A)(l - v 2 



(64) 



and we observe that we reach a point where there are 
two exactly gapless modes by tuning v — > P c , v 2 — > 1. In 
this case, the two gapless modes correspond to the char- 
acteristic excitations on an antiferromagnetic, or CDW, 
ground state, and there is no superfluid order as the 
Bloch vector is confined to the z— axis. In Tab. |TJ we 
summarize the dispersions found in the different density 
regimes in the leading order perturbation theory limit 
A = 1, and v = 0. 

In summary, we propose a conceptually simple exper- 
iment that allows to rotate the macroscopic Neel vector 
order parameter via ramping a superlattice. The mea- 
surement of superfluid and density-density correlations 
[2"2"| allows to monitor this rotation, as well as to mea- 
sure the gap of the collective pseudo-Goldstone mode, 
which is the hallmark of the proximity of the system to 
the bicritical point with enhanced symmetry. Alterna- 
tive experiments for the investigation of this proximity 
include a direct measurement of the dispersion relation 
via Bragg spectroscopy on the lattice [23], or analyzing 
the system subject to slow rotation, which also acts as a 
current defavoring SF against CDW order [2~4] . 

We further comment on the relation of our spatially 
modulated superfluid for nonvanishing v to a supersolid. 
The latter is defined as a state with simultaneously and 
spontaneously broken phase and translation symmetry. 
In our case, both symmetries are broken, but the transla- 
tion symmetry breaking is explicit and not spontaneous. 
Though the correlations are those of a supersolid, we 
would not term the state as such. 

Finally, we note that the evolution of the system from 
repulsive to attractive coupling may be viewed as a tran- 
sition from a spin 1 model (3 onsite states) to a spin 1 /2 
model. The x-y ordered phases of these two models are 
separated by the Ising transition discussed in more detail 
in the next section. 



B. Complementary exact numerical study in one 
dimension 



We now investigate how the key features of these re- 
sults manifest themselves in a ID system. This can be 
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FIG. 4. Computations of the ground state on a ID lattice 
with 60 sites for U/J = -20 using TEBD methods, (a) Cor- 
relation functions characterizing the CDW and DSF phases 
in a system with open boundary conditions ploted on a log- 
log scale as a function of distance x. Values are shown for 
N — 60 particles, where the correlation functions are almost 
identical, and N = 50 particles, where the decay of the CDW 
correlations are significantly more rapid than the DSF cor- 
relations (b) Algebraic decay exponents Kosf and Kcdw 
that are fitted to the envelope of the correlation functions for 
varying mean density n. Error bars show typical errors in the 
fitted decay. 



done by computing the ground state of the constrained 
Bose-Hubbard model using the Time Evolving Block 
Decimation (TEBD) algorithm [25 . Note that we op- 
timise our algorithm for the conserved total number of 
particles analogously to the optimisation for good 
quantum numbers in Density Matrix Renormalisation 
Group methods [27]. In Ref. [3] we already observed 
quasi off-diagonal long range order in the Dimer Super- 
fluid (DSF) correlation function {b\b\bi +x bi +x ) , together 
with exponential decay of off-diagonal elements in the 
single-particle density matrix (b.bj). This indicated the 
transition between the ASF and DSF phases in the ID 
system. 

Here we particularly investigate the situation near half- 
filling n = 1, paying attention to the interplay between 
DSF order and CDW order, characterized by the density- 
density correlation function CDW(x) = (n.jrH+x) — 
(ni)(rii +x ). In Fig. [4] we compare the DSF and CDW 
correlation functions for the ground state on a 60 site lat- 
tice with U I J = —20 and open boundary conditions. In 
Fig. |4^i we plot the correlation functions both for N = 60 
(half filling of dimers) and N = 50. At half filling the al- 
gebraic decay of these correlation functions is essentially 
the same, indeed the correlation functions are essentially 
equal, indicating coincidence of CDW and DSF orders in 
this state. Whilst reducing the total number of particles 
on the lattice to N = 50 does not significantly change 
the DSF, the density-density correlation function decays 
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FIG. 5. Numerical validation of symmetry enhancement 
SO(2) — > SO(3). The plot shows the exponents Kdsf and 
-Kclw describing the algebraic decay of DSF and CDW or- 
der in the ground state as a function of the ratio U / J at unit 
filling n — 1. For sufficiently large interactions, the coinci- 
dence of the decay exponents signals the degeneracy of the 
two kinds of order. These calculations were performed with 
TEBD methods for 60 particles on 60 sites. Open bound- 
ary conditions were used, but the decay exponents were fitted 
within the central 30 sites. The fitting errors are similar to 
those in Fig. [4] (b). Number of Schmidt coefficients retained 
in TEBD calculations x = 200. 



much more rapidly in the ground state, in addition to 
large superimposed oscillations. This relative sensitivity 
of the correlation functions is characterised in Fig. |4]d, 
where we show the result of fitting an algebraic decay 
x * to the envelope of each of the correlation functions. 
Again, we see that the decay of CDW and DSF correla- 
tions is identical within fitting errors at unit filling, but 
the CDW is very sensitive to deviations from unit filling, 
and it is dominated away from n = 1 by the DSF. 

In Fig. [5] we study the approach of the bicritical point 
at fixed half filling n = 1 as a function of the ratio of 
hopping and interaction J/U. The plot clearly shows the 
symmetry enhancement from the conventional U(l) ~ 
50(2) to an SO(3) symmetry: The two decay exponents 
describing DSF and CDW order approach each other for 
sufficiently strong attractive onsite interaction U, thus 
indicating the degeneracy of the two different kinds of 
order. 

In an experiment it would be difficult to produce a 
setup with an exactly commensurate number of parti- 
cles and lattice sites. One way to observe emergence of 
the CDW order, though would be to prepare the sys- 
tem in a harmonic trapping potential, where the density 
would vary across the trap. In Fig. [6] we investigate the 
ground state for 30 particles on 60 lattice sites in the 
presence of such an external harmonic trap. In Figs. [6^,, b 
we show the correlation functions for CDW and DSF or- 
der as they vary across the trap. We note that the DSF 
order is significant throughout the occupied region. We 
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FIG. 6. Computations of the ground state with 30 particles 
on 60 lattice sites with U/J = —20, in the presence of a har- 
monic trap with on-site potential V(x). (a) Shaded plot of 
the CDW correlation function (n x n y ) — (n x ){n y ) (with inter- 
polated colours, and diagonal elements not shown), showing 
substantial order near x = y — 20 and x = y = 40, where the 
mean filling factor n ~ 1 regions. V^a;) = Vt T (x — 30. 5) 2 , 
Kr = 1.33 x 10~ 3 . (b) Shaded plot of the DSF correla- 
tion function DSF(x — y) = (b x b x b y b y ) (with interpolated 
colours), showing substantial order across the occupied re- 
gion of the lattice, for the same trap parameters as part 
a. (c) Density n(x) for the same parameters as parts a,b. 
(d) Plots of the DSF correlation function DSF(x) for a trap 
V(x) = Vt x (x - 30) 2 , Kr = 0.7/30 2 « 7.78 x 10" 4 , with and 
without an additional superlattice potential )*. 



have checked in addition that across this region, the off- 
diagonal elements decay algebraically as a function of 
distance. On the other hand, the CDW correlations are 
most significant in regions near unit filling. For the trap 
parameters chosen here, this occurs near sites 20 and 40, 
as shown in Fig. [6J3, where we also see significant oscilla- 
tions in the density, which are also characteristic of the 
appearance of CDW order. In Fig. [6]i we then investi- 
gate how the order can be manipulated by the addition 
of a weak superlattice. We see that the addition of an al- 
ternating potential on the order 0.01J is sufficient to sig- 
nificantly increase the algebraic decay exponent for DSF 
order. Because the system size is small, it was difficult to 
obtain reliable results for the algebraic decay exponent of 
CDW order, but our calculations indicate that applying 
such a superlattice can indeed be used to select the dom- 
inant order for a system in the presence of a harmonic 
trap. 

Using t-DMRG methods we can also investigated pos- 
sible time-dependent preparation of the continuous su- 
persolid beginning from a Mott insulating state in the 
presence of a superlattice, analogously to the studies per- 
formed in Ref. 3J. Beginning in an insulating state with 
two atoms in the lowest wells of a period two superlat- 
tice, it is possible to prepare a state with n = 1 and 
U/J = —20 in a timescale of the order of 100 J -1 , with 
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good fidelity of the DSF correlation functions, provided 
that a sufficiently strong constraint can be imposed so 
that no loss events occur on the timescale of the ramp. 



be lead based on a simple quadratic spin wave theory. 



A. Low Energy Derivative Expansion 



V. LONG WAVELENGTH LIMIT: NATURE OF 
THE PHASE TRANSITION 

Even at low energies, non-linearities in the effective 
action may in principle have an impact on the physi- 
cal observables, such as the nature of the phase tran- 
sition. Such a scenario is known as Coleman- Weinberg 
phenomenon [T3]: Two near gapless degrees of freedom 
are coupled to each other, in a way that a phase tran- 
sition which one of them undergoes is driven first order 
due to the long wavelength fluctuations of the other: the 
first order transition is radiatively induced. 

In our problem, indeed we face competing low energy 
degrees of freedom at the ASF-DSF transition: First, 
there is the gapless Goldstone mode present in the dimer 
condensate, which does not undergo qualitative changes 
at the ASF-DSF transition point. Second, at the Ising 
type transition one expects a Z 2 degree of freedom to 
emerge in the low energy sector for the atom degrees 
of freedom. A possible coupling between those low en- 
ergy degrees of freedom may or may not give rise to a 
Coleman- Weinberg mechanism. 

Here we study this question by means of a systematic 
derivative expansion of the effective action. At "low" 
densities n rs 0, 2 we identify a first order transition 
in line with known results for continuum bosonic Fes- 
hbach models at low density [BHH]. Such a reproduction 
of the continuum results must be generally expected in 
low density lattice systems. This situation is seen to be 
rather generic in nonrelativistic systems [28-30 . But, 
intriguingly, there is a lattice based decoupling mecha- 
nism which guarantees the existence of a second order 
transition, and thus a true quantum critical point, in the 
vicinity of n = 1. Thus, we identify a true Ising quantum 
critical point in our system, connecting the two ordered 
ASF and DSF phases. 

Note that the scenario crucially hinges on the control 
over a coupling of the two near gapless modes close to 
the transition. It is evident that the discussion cannot 



Our strategy is as follows: We will approach the phase 
transition from the DSF side, where there is not yet an 
atomic condensate, and tune the atomic mass parame- 
ters to criticality from there. For this purpose, we draw 
the low energy, continuum limit of the effective action 
corresponding to Eq. ( 30 ) . We then identify the relevant 



low energy fluctuations and integrate out the massive de- 
grees of freedom. We arrive at an action that describes 
the dimer Goldstone physics, the Ising degree of freedom 
as well as a cubic coupling of Goldstone mode to Ising 
density. The derivation is similar to the one presented 
in [7] in the continuum, differs however in the crucial 
aspect that the model discussed there features already 
microscopic propagating dimer degrees of freedom. Here 
we show how such terms are generated via successive in- 
tegration of the massive degrees of freedom. 

At low energies, the action corresponding to the Hamil- 
tonian (30) encounters two immediate simplifications: 
First, we consider the constraint X+ = 1 — hu — h%^: 
The density operators are less relevant than the number 
1 at low energy. Consequently we replace Xi — > 1. By 
this replacement, we effectively drop the local constraint 
for the atoms and dimers. Physically, this is justified 
from the fact that infrared fluctuations with wavelengths 
much larger than the lattice spacing do not resolve single 
sites - as stated above, while symmetries provide scale 
independent restrictions on the form of the effective ac- 
tion, the relevance of the constraint principle depends on 
scale. Second, we draw the continuum limit. Our origi- 
nal Hamiltonian (30) often contains bilocal terms. In the 



quadratic sector, the resulting spatial derivative terms 
are kept: they describe spatial propagation and may be 
leading in the infrared for zero mass terms encountered 
close to the phase transition. However, in the interac- 
tion terms we drop the gradient couplings if they appear 
in combination with a local one, which in comparison is 
always more relevant in the sense of the renormalization 
group. Finally, we drop the local quartic terms coupling 
atoms with dimers, which are subleading in comparison 
with the local cubic ones. The corresponding action reads 



S = Si [b\ , h] + S 2 [b\ , b 2 ] + S int [b\ , h , a, tt] , (65) 
5i[6i] = f b\[d T -n- (-2/1 + U)s 2 ) - J(l + |s| 2 )(z + A)]6i - V^Jcs^z + A)6j + c.c), 

J X 

S2M = I b\[d T + (-2/1 + U)(c 2 - s 2 )]b 2 , 



S Int [61, 6a] = / Jz[3cs {b\+b 2 )n 1 - V2(c 2 - s 2 )(b\b\ + c.c.)] . 

J x 



(J = J d,Td d x,x — (t,x),A the Laplace operator. We omit the (r, x) dependence of the field for brevity.) Here 



18 



and in the following we have chosen s real without loss 
of generality. 

In the next step we identify the relevant phase fluc- 
tuations. The terms in the action (65 1 are seen to be in 



two classes: The first one is made up of field combinations 
which transform according to 17(1) x 17(1) (the first phase 
rotation acts on bi and the second on & 2 ); i-e. they do not 
lock the phases. In contrast, the cubic interaction terms 
in the last line of Eq. ( [65] ) lock the phases such that the 
residual symmetry is a single 17(1). (Such a mechanism 
breaking 17(1) x 17(1) — > 17(1) is a consistency check for 
our theory, which emerges from a constrained version of 
the Bose-Hubbard model, in turn only possessing a sin- 
gle 17(1) symmetry.) The dominant temporal and spatial 
phase fluctuations thus originate from the vicinity of the 
phase constraint emerging from the phase locking of the 
atomic to the dimer phase, 62 (x) = 20i(x). To bring 



out the physics of these fluctuations, it is convenient to 
perform a local gauge transformation on the b\ field such 
as to absorb the 9 2 fluctuations [7]. Here we work in 
cartesian coordinates for the fluctuating fields, and con- 
sequently the gauge transformation is realized linearly. 
The gapless phase fluctuations of the dimer field are rep- 
resented by its imaginary part, b^ix) — (a(x)+in(x))/\/2 
(cf. Eq. (15)). To absorb the phase fluctuations into b\, 



we introduce dressed fields according to 

b x (x) -> bx(x) = h(x)(l - iKir(x)), (66) 
k = (c 2 - s 2 )/(2V2cs). 

Now the gauge transformed action can be calculated. In 
this expression, we only keep leading terms which are 
affected at linear order in the infinitesimal rotation. The 
result is 



Si[b\M 



2 l(b\M) 



d T + m\ — J(l + s 2 )A 
-2y/2Jcs(z + A) 



-2V2Jcs(z + A 



-Br 




(67) 



Sint \y\ 1 b\ > 7r] = / V2iK,d T irb\bi + Jza 



3V2cs b\bx - (c 2 - s 2 )(6ife! + b\b\) 



with ml w \U\/2 - Jz(l + s 2 ), using -/z w \U\/2 - as 
we are only interested in the low energy limit, the precise 
value of the couplings is unimportant, and we will work 
with the mean field values (which are, however, expected 
to be rather accurate except for the small density regime 
n ~ 0, cf. Sec. III). As a preparation for the elimination 
of the massive modes, we further introduce hermitian 
fields for the single particle excitations b\ = tp + hp, such 
that the action reads 



Si [<p,ip] 



-\d T m 2 _ — ^ 2 A 




5 int [^,V,cr,7r]=/ i^d T ir(<p 2 + i/> 2 ) + o-(A_^ 2 + A+^ 2 ), 



where £± = J(c ± -\/2s) 2 , m\ 



(68) 



Jz(3cs/V2 ± (c 2 - s 2 )), and again we keep only leading 
terms. As appropriate for the phase mode, the field tt 
interacts with the atomic fields ip and ip only through 
its time derivative, while the field a interacts directly. 
Note that m\ > m_, and upon approaching the phase 
transition, m + hits zero prior to m?_ 51 [7J . Indeed the 
condition ml = coincides with Eq. ( 45 ) if we also use 
the mean field equation of state n = 2s l . For vanishing 
m+, we then find m 2 _ A\/2csJz = Ay/n(l — n/2)Jz 
within the mean field approximation for the high energy 
physics. Hence, the field i/j (the imaginary part of the 
atomic field b\ ) remains massive for any density < n < 
2 at the transition, and we may safely integrate it out 



perturbatively at the one-loop level, while the remaining 
degree of freedom ip becomes soft and plays the role of 
an Ising field. The resulting effective action for the fields 
ip, tt, and a reads 

S[<p,TT,a]=J |^o-(A/ 2 -C 2 A)cr + icra r 7r + i-^a r 7r^ 2 
+((d T Tr) 2 - £ 2 7rA7r 

+ l -ip{m 2 + -Z v d 2 T -e + ^}, (69) 

where M 2 - X 2 _/m 2 _, C ~ K 2 /m 2 _, £ 2 - £ 2 - A 2 ^ 2 /mi, 
and Z v ~ mZ 2 - (Note that in the limit |17| ^> Jz both 
fields ip and ip are massive and, after integrating them out 
perturbatively, we get Eq. (46) for the effective dimer 
Hamiltonian of the 62 field.) The field a now becomes 
massive and can be integrated out as well. The final 
effective action for the fields ip and tt is 

s eS [<p,TT] = j x {y(-z v d 2 - e+ A + m 2 + V + 

+ ^Zd 2 -eA)n + i^d T np 2 }, (70) 

with Z ~ M~ 2 and A ~ A^/M 2 . This action describes a 
coupled theory for the Goldstone mode w and the Ising 
mode ip. Note that here we also keep a fourth order Ising 
coupling. Its presence being rooted in the tree-level a 
exchange, this coupling is positive. Thus, the low energy 
theory contains an Ising part, i.e. a real field with quartic 
potential which exhibits Z 2 symmetry breaking when m 2 ^ 
turns negative. If this part of the action were isolated, 
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the transition would be in the Ising universality class, 
and therefore of second order. In the presence of the 
Goldstone-Ising coupling, more care needs to be taken: 
In general, a coupling of two (near) gapless real bosonic 
degrees of freedom can lead to a fluctuation induced first 
order phase transition, known as the Coleman- Weinberg 
phenomenon [15] . The Ising self- interaction A and the 
Ising- Goldstone coupling k can be compared via naive 
power counting the canonical dimension of A is 3 — d, 
and that of k is (3 — d)/2. Thus, in any dimension the 
corresponding terms have the same degree of relevance 
and therefore compete with each other. 



The form of the action ( 70 ) coincides with the one ob- 
tained in [7] from the continuum Feshbach model. The 
renormalization group analysis of the action ( 70 1 for 
nonzero k has been performed in d — 3 by Frey and 
Balents [SB] at T = 0, and extended to nonzero tempera- 
ture by Lee and Lee [2S], revealing a Coleman- Weinberg 
phenomenon. Thus, for a generic k^O the phase transi- 
tion will be driven first order. This scenario is realized in 
the low density limits n rs 0, 2, where our conclusion thus 
matches the expectations from the continuum, which was 
anticipated in [HJ [8] and discussed in detail in [7] . 

However, the lattice offers the possibility to penetrate 
the regime where n ss 1. Here, an intriguing situation 
appears: There exists a point in the phase diagram at 
which the coefficient of the cubic terms vanishes exactly, 
which happens du e to the zero crossing of the coupling 



k. From Eq. (66) we have k ~ c — s . Working with 
the mean field equation of state n — 2s 2 , one concludes 
that this takes place at n = 1. In reality, renormaliza- 
tion effects will add contributions to the naive value of 
k. Furthermore, inspection of the full equation of state 



( 25 1 suggests further shifts from the naive expectation, 
but we have seen in Sec. IIII Bl that close to n = 1 these 
are small. Thus, we expect the decoupling point to be 
located in the close vicinity of the commensurate point 
n = 1. We provide further evidence for this expectation 
from a symmetry argument in the next section. 



B. Symmetry argument for the Ising quantum 
critical point 

The decoupling of Goldstone and Ising mode at a spe- 
cial point in the phase diagram can also be obtained from 
a symmetry argument. Being based on a combination of 
the phase locking symmetry between the degrees of free- 
dom bi , b 2 and a temporally local gauge invariance, it 
complements the above explicit derivative expansion and 
sheds more light on the origin of the decoupling of Ising 
and Goldstone physics. 



4 The power counting applied here is based on the effective rela- 
tivistic Ising and Goldstone low energy actions with dynamical 
exponent z = 1, and not the original nonrclativistic theory 



For this purpose, let us first discuss the temporally lo- 
cal gauge invariance of the Bose-Hubbard Hamiltonian 
[31] in the presence of an infinite three-body repulsion, 
which is equivalent to the constrained model under con- 
sideration here. This adds a local term to the standard 
Bose-Hubbard Hamiltonian, 

H= lim [H BH + 73 V;n i (n i -l)(n i -2)]. (71) 



The temporally local gauge invariance results simply 
from the fact that the Hamiltonian is not explicitly time 
dependent (while it is spatially non-local, such that a 
spatially local gauge invariance does not exist). Conse- 
quently, the constrained Bose-Hubbard action must take 
the form 

Sc = J dTY^(al(d T -ti)ai + H[a\a]) (72) 



such that the temporally local gauge invariance is ex- 
pressed as an invariance under 

a, — » cxpiA(t)ai, /i — > id T X(t). (73) 

Since our construction must conserve this property, we 
also require this invariance for the theory defined with 
(30 1. On the level of the effective action and in Fourier 



space, this invariance translates into the Ward identity 
for the effective action 



d 8 2 T 
~d»6b\ /2 (q)5b 1/2 (q) 



b 1/2 =0;q=0 



(74) 



d 



S 2 T 



9 M Sb\ /2 (q)6b 1/2 (q) 



bi/2=0;<J=0 



i.e. the coefficient of the linear time derivative must 
equal the derivative with respect to the chemical poten- 
tial. Therefore, in a derivative expansion of the effective 
action, which is appropriate at low energies, we have: 



+ m'i + ...]&! +^1 2 + b'i) 



TJb\[z 1 d T + Vl d' i T 
+b\[z 2 d T + y 2 d 2 + m 2 2 + ...]b 2 + h(b\b\ + b 2 b\ z ) + . 



(75) 



The presence of a condensate for b 2 , 9 ^ 0, generates off- 
diagonal terms in the b\ inverse propagator, i.e. £ =/= 0. 
Here we restrict to the spatially local part of the effec- 
tive action, since this is the sector where the coupling 
of Ising to Goldstone mode emerges. The Ward identity 



(74} implies z 1/2 = ~dm 2 /2 /dii =: gxji 



Furthermore, using solely the global gauge invariance, 
we can make the connection between g 2 and g\. Indeed, 
we have a phase locking in the i^( 21 )i^( 10 ) t _i_ h.c. term. 
As a consequence of these terms, the phases of b\ and b 2 
cannot transform independently, and we have 



bi i expiA&i i, bn —¥ exp 2i\b 2 



(76) 
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leading to the additional Ward identity 
_ d 5 2 T d S 2 T 



dn Sb\{q)Sbi(q) b 1/2 =o, q =a dn 5b 2 (q)Sb 2 {q) 



6i/2=0;?=0 

(77) 



or .92 = 2gi- In sum, we have the following relations: 



Z-2 



92 = 2^1 = 20i. 



(78) 



Next we discuss properties of the "compressibility" 
coupling g2{n) = —dm\/dy,\ n , which fixes how strongly 
the bound state excitation couples to the chemical po- 
tential. In the limits n ~ 0, 2 we can compute it exactly 
from the solution of the corresponding two-body prob- 
lems Eqs. (33). At n = 0, we find g 2 > 0, while at n = 2 



we obtain g2 < 0. These opposite signs can be expected, 
as at n = the excitations are well-defined dimers, while 
at n — 2 we face well defined di-holes. If we do not re- 
define the chemical potential, then adding a di-hole is 
energetically equivalent to delete a dimer. Under the 
mild assumption that the compressibility is a continuous 
monotonic function of n (our description is tailored to 
describe the DSF phase including the phase border, and 
therein we do not expect additional phase transitions), 
then g 2 (n) must have a unique zero crossing. We note 
that we should use the above derivative prescription as 
an operational definition of g 2 ; in principle, there could 
be a /^-independent constant adding to the full mass or 
gap term of b 2 . For b\, such a situation takes actually 
place and we have an additional mass or gap term U. 

As a consequence of Eq. ([78]), a zero crossing of g 2 
also implies a zero crossing of the coefficients z 2 , g\ , z\ . 
Thus, the leading frequency dependence is not linear, but 
quadratic, and the analogous statement is valid in the 
time domain, where the leading behavior is a quadrati- 
cally appearing time derivative. 

With this result, we now discuss the possible form of 
the coupling of the Ising to the Goldstone mode. As 
above, we decompose linearly into massive and phase 
mode, and absorb the phase fluctuations into dressed b\ 
fields, bx ->■ bx = &i(l - ikn/2 3 / 2 ), k = h/L Indeed the 
low energy effective action can only depend on derivative 
couplings associated to the phase mode tt, due to the 
global U(l) invariance under transformations ir — > n + A. 
The transformaton cancels the cubic term in Eq. ( 75 ) 



associated to phase fluctuations, while the contribution 
associated to the real part a can be dropped at low en- 
ergies since the amplitude is massive. At the same time, 
the b± part in the dressed frame now reads 



b\[z l d T +y l d 2 + m{ + ..]b l 



(79) 



-izikd T Trb\bi + iyikd 2 irb\~bi 



Thus, for g 2 = 0, Eq. (78) also implies that the cu- 



bic derivative coupling z\k with canonical dimension 
(3 — d)/2 vanishes. The leading term is a cubic cou- 
pling with quadratic time derivative. This coupling has 
canonical dimension (1 — cZ)/2, and thus is irrelevant near 



a Gaussian fixed point for d > 1. Similarly, a potential 
U(l) symmetric coupling term g' f(d T ir) 2 (j) 2 has canon- 
ical dimension 1 — d. Both therefore do not lead to a 
Coleman- Weinberg phenomenon. In consequence, Gold- 
stone and Ising physics effectively decouple at low ener- 
gies, giving rise to a second order Ising transition. 

We summarize our result. Based on the zero crossing 
of g 2l phase locking and temporally local gauge invari- 
ance we find: 

(i) At the zero crossing point, the nonrelativistic time 
derivative terms vanish. In the sense of a derivative ex- 
pansion, the next relevant term is d 2 , in which case the 
theory acquires a relativistic space-time isotropy in a d+1 
dimensional space-time. This is physically sound, as this 
point has a special kind of (di-)particle-hole symmetry, in 
that the hybrid excitation consists of a superposition of 
"dimers" and "di-holes" to equal parts. However, we note 
the absence of a particle-hole symmetry in the conven- 
tional sense - such a situation only occurs in the pertur- 
bative limit J/\U\ —> 0, as discussed in Sec. IV Beyond 
the leading order perturbation theory, this symmetry is 
broken. One manifestation of the absence of this sym- 
metry is the asymmetry of the critical line in the phase 
diagram, cf. Fig. [2j 

(ii) The cubic coupling of Goldstone to Ising mode also 
vanishes at this point. Only terms which are irrelevant 
in d > 1 dimensions then can couple these modes. As a 
consequence, the Coleman- Weinberg mechanism is sup- 
pressed. 

We observe that the constraint influences the physics 
even at very long wavelengths: It is responsible for the 
existence of a maximum filling, in turn leading to the 
existence of a zero crossing of the dimer compressibility, 
in turn responsible for the existence of the Ising quantum 
critical point. 

In conclusion, close to the "particle-hole symmetric" 
point at n = 1, there is a d + 1 dimensional Ising quan- 
tum critical point. Examples of physical realizations of 
Ising quantum critical points in nature are actually rare. 
Several systems exhibit Ising type phase transitions with 
discrete symmetry breaking, like the ASF-DSF transition 
in the continuum Feshbach model [5] and or a transi- 
tion between superconductors with different pairing sym- 
metries [3D], but in these cases in the long wavelength 
limit a Coleman- Weinberg phenomenon takes place. A 
cubic coupling of the Goldstone mode with linear time 
derivative to the Ising density is actually quite generic in 
nonrelativisic systems, where the Ising mode emerges as 
an effective degree of freedom describing the transition 
from one ordered phase to the other. Here we have iden- 
tified a mechanism that suppresses this coupling. One 
of the few other examples for Ising quantum criticality 
is possibly provided by the model magnet LiHoF4 [32], 
though the issue is debatable due to the long range in- 
teractions in the material, preventing an exact mapping 
to the Ising model. 

The fact that qualitative aspects of the critical behav- 
ior are changed in the vicinity of the particle-hole sym- 
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metric point n = 1 bears some resemblance to the physics 
at the tip of the Mott lobe in the repulsive Bose- Hubbard 
model. There, the behavior changes from the nonrela- 
tivistic 0(2) (or XY) universality class with dynamical 
exponent z = 2 to the relativistic 0(2) model with z = 1 
(331. 



C. Estimate of the Correlation Length 

To get an impression of the perspective to observe Ising 
quantum criticality in this system experimentally we es- 
timate the correlation length. This quantity is accessible 
with current experimental technology |34| . and has been 
measured in continuum Bose gases to characterize critical 
behavior. 

The Coleman- Weinberg phenomenon manifests itself 
in the presence of "runaway" trajectories on the RG 
flow diagram. We therefore can estimate the correlation 
length at the first order phase transition as a scale Z*, 
at which the runaway trajectory with the corresponding 
initial conditions hits the boundary of the stability re- 
gion of the system [IT]- The instability is characterized 
by the quartic Ising coupling A turning negative, i.e. the 
condition A (Z*) = 0. 



The scaling properties of the action (701 are deter- 
mined by three parameters: V = (£/£ + )y/Z v /Z, U = 
and K = K 2 i 2 /SXZ^/Z~^ with the corre- 
sponding RG equations derived in Ref. [28]. The quan- 
tity V scales to zero, therefore we can put V — from 
the very beginning. Then the RG equations for the re- 
maining constants K and U read 



1 dK 1 TT 5 ^ 

K~dT =£ ~ 4 2 



1 dU 
U~dJ 



-U -6K -24 



U : 



(80) 
(81) 



where e = 3 — d. To solve these equations, we first in- 
troduce new functions k = K exp(-el), u = U exp(— eZ), 
and a new variable x = exp(eZ). The equations then have 
the following form, 



dk 
dx 
du 
dx 



1 , 5 ,2 

-uk + -At 
4 2 



6uk + 24k' 



(82) 
(83) 



Writing u = kf(k) and, therefore du/dk = f + kf, we 
obtain 

df _ du/dx _ f 2 + 28/ + 192 _ (/+ 12)(/ + 16) 
k— — ,. , . ,/ — 



dk dk/dx 



2/ + 20 



2(/ + 10) 



This equation can easily be solved with the result 



k 




-i 6 y 




f 12\ 


k 


V/o- 


f 16/ 




-12 J 



(84) 



(85) 



where /o = Uo/ko is the initial value for the function / 
when k = kg. 



It follows from Eq. (831 that 



£ | = -¥[ S / 2 + 48/ + 192] 



e Tx [k fm 



dx 



f dk J.U 

d f 



and, after using Eq. ( 84 1 , we obtain 
df 



dx 



-(/ + 16)(/ + 12) 
i fco(/o + 12)(^^ 



(Jo + 16) 3 
The solution of this equation reads 

x — 1 _ exp(eZ) — 1 8 1 



3fc f + 12 



1 - 



fo + 16 
/ + 16 
(86) 

and, together with Eq. (85), provides a general solution 
of the RG equations (82) and (83) and, therefore (80) 



and (81) 



The above solution allows us to find the scale Z*, at 
which the RG flow reaches the border of stability, U(l*) — 
0. In 3D we obtain (after taking the limit e = 3 — d — » 0) 



8 



3fc (/o + 12) 



1- 1 



fo 
16 



with fo = uo/ko. As a result, close to the Ising critical 
point, fc 0, we get 



(1-n) 



This result indicates a rather broad critical domain in 
density around the true Ising critical point, in which the 
correlation length extends over the whole system. Such 
extended quasi-critical behavior can be expected for a 
fluctuation induced first order transition, which results 
exclusively from the competition of very long wavelength 
degrees of freedom, and therefore should be weak. For 
example, already at filling n = 1/4, 2 — 1/4 the correla- 
tion length is on the order of 15 lattice sites, and greatly 
exceeds the typical size of an optical lattice of 20 to 100 
sites at filling 1/2, 3/2 already by a factor of 10. We con- 
clude that the Ising quantum critical behavior should be 
experimentally observable in our system. 

Finally, we emphasize that the discussion presented in 
this section crucially hinges on the fact that our field 
theoretic setup allows to fully assess the effects of inter- 
actions, i.e. nonlincarities in the effective action. Here 
we have shown that these interaction effects persist even 
down to arbitrarily long wavelengths. Obviously, such a 
scenario is not captured in a simple quadratic spin wave 
theory with a priori decoupled atomic and dimer excita- 
tions. 
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VI. CONCLUSION 

In this paper, we have performed a detailed analyti- 
cal investigation of the phase diagram of the attractive 
lattice Bose gas with a 3-body hardcore constraint. For 
this purpose, we make use of a method presented in [11] 
which allows to exactly map the constrained model to a 
theory for two unconstrained bosonic degrees of freedom 
with conventional polynomial interactions. Within this 
framework, we particularly focus on effects tied to inter- 
actions, which cannot be addressed within a mean field 
plus spin wave approach. While our analysis confirms the 
rough features of the phase diagram obtained from a sim- 
ple mean field approach - the presence of an Ising-type 
phase transition from an atomic to a dimer superfluid, 
numerous interaction driven effects are identified. These 
arise on various length scales, ranging from the fluctua- 
tion induced formation of the dimer (or di-hole) bound 
state on top of the vacua at n — and n = 2 on the 
microscopic level over a an understanding of the beyond 
mean field effects causing nonuniversal shifts in the phase 



boundary and giving rise to the proximity of the system 
to a bicritical point with enhanced 50(3) symmetry in 
strong coupling, down to the assessment of the true na- 
ture of the phase transition at very long wavelength. This 
underpins the fact that short and long range correlations 
can then be treated within a unified formalism. 
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